Lagrange mesh and exact diagonalization for numerical study of semiconductor 
quantum dot systems with apphcation in singlet-triplet qubits 
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We present a highly flexible computational scheme for studying correlated electrons confined by an 
arbitrary external potential in two-dimensional semiconductor quantum dots. The method starts by 
a Lagrange mesh calculation for the single-particle states, followed by the calculation of the Coulomb 
interaction matrix elements between these, and combining both in the exact diagonalization of the 
many-body Hamiltonian. We apply the method in simulation of double quantum dot singlet-triplet 
qubits. We simulate the full quantum control and dynamics of one singlet-triplet qubit. We also 
use our method to provide an exact diagonalization based first-principles model for studying two 
singlet-triplet qubits and their capacitative coupling via the long-distance Coulomb interaction. 

PACS numbers: 73.22.-f,81.07.Ta 
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I. INTRODUCTION 

The development of experimental methods has enabled 
the fabrication of "artificial atoms" with a controlled 
number of electrons, ranging from a few to a few hun- 
dred, confined in a tunable external potential inside a 
semiconductor.^"^ These quantum dots (QD's) have been 
proposed as a possible realization for the qubit of a quan- 
tum computer^'^. 

A framework for using two-electron spin eigenstates 
as qubits was proposed by Levy in 2002^. The two- 
electron double quantum dot (DQD) spin states have nat- 
ural protection against the decoherence by the hyperfine 
interaction and allow a scalable architecture for quan- 
tum computation^. The universal set of quantum gates 
for two spin singlet-triplet DQD qubits has been demon- 
strated experimentally. These gates include one qubit 
rotations generated by the exchange interaction^ and sta- 
bilized hyperfine magnetic field gradients^, and two qubit 
operations using long distance capacitative coupling by 
the Coulomb^'^ interaction. 

In creating the inter-qubit gates and operations, quan- 
tum entanglement is essential^ ^. The aforementioned 
capacitative coupling is one possible method to create 
entangled states and implement two-qubit operations 
in singlet-triplet qubits, the other possibility being ex- 
change based methods®'^^. In the capacitative dipole- 
dipole coupling, the entanglement is achieved by differ- 
ing charge densities in the singlet and triplet states that 
result in different Coulomb repulsion between the qubits. 
This conditioning can be used to create entangled states 
and to implement the two-qubit gates required for uni- 
versal quantum computing'''^^'^'^"^^. 

Although other methods, like the variational quantum 
Monte Carlo^® and the density functional theory^^, have 
shown to give reasonably accurate results, exact diago- 
nalization is still the most reliable technique for small 
particle numbers. In this paper, we use the Lagrange 
mesh method^* and exact diagonalization to simulate 
one and two singlet-triplet qubit systems. The one- 



particle eigenstates are computed using the Lagrange 
mesh method and are then used to create the many-body 
basis in the calculations. The use of the Lagrange mesh 
method as a source of single-particle states allows us to 
study very flexibly various forms of the confinement po- 
tential. For example, the relatively complex two-DQD 
case can be handled with ease using this method. Still, 
the one-body basis needed for good accuracy is much 
more compact in the Lagrange mesh than, e.g., in the 
finite-difference formulation. 

This paper is organized as follows. In Section II, the 
theoretical model used in our computations is briefly 
discussed. In Section III, we introduce the Lagrange 
mesh method for many-body problems in quantum dots. 
The computational results are shown in Section IV. The 
single-particle states and their convergence are discussed 
in Section IV A. In Section IV B, we model the full quan- 
tum control and dynamics of a singlet-triplet qubit. In 
Section IV C, we use the Lagrange mesh method to cre- 
ate a realistic first-principles ED model for studying the 
interplay and entanglement of two singlet-triplet qubits. 



II. MODEL 

We model a lateral GaAs quantum dot system with 
the two-dimensional Hamiltonian 



N r 



H 



E 



2m* 



+ V{r,) + Vz{rj) 



j<k 



'jk 



where Vz{rj) — g* fiBB{rj) ■ Sj is the Zeeman term 
with the effective GaAs g-factor g* = —0.44. A is 
the magnetic vector potential, and to* « 0.067 mg and 
e ~ 12.7 eo are the effective electron mass and permittiv- 
ity in GaAs, respectively. In numerical work, it is con- 
venient to switch into effective atomic units by setting 
TO* = e = h = l/Ane = 1. In these units, energy is given 
by Ha* ss 11.30 meV and length in Oq w 10.03 nm. 
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In our computations, the external potential V{r) for 
quantum dot systems consists of several parabolic wells. 
A confinement potential of n parabolic wells can be writ- 
ten as 



V{r) = -m*ojl min {\r - r^^}, 



l<j<n 



(1) 



where {'"i}i<j<ri are the locations of the minima of the 
parabolic wells, and ljq is the confinement strength. The 
kinks caused by the min-function are smoothed in the 
potential. 



III. METHOD 

The Lagrange mesh method^^ is a very efficient 
method for solving the Schrodinger equation, and it has 
the simplicity of a finite-difference mesh calculation, since 
no integrations need to be performed. It also does not 
suffer from the same limitations regarding to the con- 
finement potential as for example using the analytical 
Fock-Darwin basis. 

We will use this technique to solve the eigenstates of 
the one particle Hamiltonian, 



H = + V{r) 
2m* ^ ^ 



(2) 



omitting the magnetic vector potential here. The eigen- 
states are then used as a basis for the many-body calcu- 
lation done by the exact diagonalization technique, and 
the Zeeman term and Coulomb interaction are included 
in it. 



A. One-particle problem 

A set of N Lagrange hmctions L/. defined over an in- 
terval (a, b) is associated with N mesh points Xk & {a, b) 
and a corresponding Gauss quadrature 



nb N 

/ dxf{x) K^\kf{xk) . 

k=i 



(3) 



The Lagrange functions are infinitely differentiable real 
functions, which are orthonormal, 



/ AxLi{x)Lj{x) = Sij , 

J a 

and satisfy the Lagrange conditions 



(4) 



(5) 



From the conditions of Eqs. (4) and (5), it follows that 
the Gauss quadrature is exact for any product of two 
Lagrange functions: 

j dxLi{x)Lj{x) = 5ij = ^ \kLi{xk)Lj{xk) ■ 

k=i 



Many different Lagrange meshes, mostly based on or- 
thogonal polynomials or trigonometric functions, have 
been proposed^^ both for finite intervals and the infinite 
intervals (0,oo) and (—00,00). The meshes can also be 
modified to distribute the mesh points optimally for a 
particular system.^" 

One of the most simple Lagrange meshes is the sine 
mesh.^^ It is defined over the interval (—00, 00), but de- 
signed to treat fairly well localized wave functions. The 
mesh points distributed uniformly around the origin are 

xa = a, ae{-^,-^ + l,...,^} , (6) 

and all the weights in the Gauss quadrature are = 1. 
The Lagrange-sinc functions are 



La{x) = sinc(a; — a) = 



sin[7r(x — a)] 
7r(x — a) 



The matrix elements of the derivatives dx and be- 
tween two sine functions can be calculated analytically, 
resulting in 



and 



/ dxLa' 

J — 00 


{x)dxLa{x) 
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The potential energy matrix elements can be calculated 

analytically for some potentials, but it turns out that for 
the smooth potentials, these can be accurately approxi- 
mated using the Gauss quadrature of Eq. (3) as 

Va'a « V{a)Sa'a ■ 

Strictly speaking, this approximation breaks the varia- 
tional principle. The validity of the Gauss quadrature 
approximation is discussed in the Appendix. 

Generalized to two dimensions and an area 
^ ¥)' Lagrange-sinc functions are 

given by 

La{r) = f sine [f - x^J] sine [f (2/ - , (7) 
where the N x N mesh points are scaled to 



' a — N 



a , a^, tty e {- 



Af-l N-1 
2 ' 2 



JV-1 



1 , . . . , 2 



having grid spacing h = L/N and weights Xg. = h^- The 
matrix elements of the Hamiltonian (2) between the basis 
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functions in Eq. (7), in effective atomic units, are 



' 5^ + V{ra) , a'x = ax,a'y = ay 



H, 







a'x ¥= a,x,a'y = ay 



, a'^ ax,a'y ^ ay , 



where V{r) is the external potential. Diagonalization 
of this Hamiltonian matrix gives the one-particle eigen- 
states and -energies. The accuracy of the results obtained 
can be tested by varying the number of mesh functions 
N and the side length of the simulation square L. 



where the expansion coefficients a multiply the inter- 
action matrix elements Vabcd between the sine basis 
functions. To calculate these, we start with the two- 
dimensional Fourier transform of l/ri2, namely 



J" 



1 

ri2 



(fe) = / d 



ri2- 



ri2 

27r 



Jo Jo 



2 cos(0-t-7r) 



where 6 is the angle between and k. Using the Jacobi- 
Anger identity of Bessel functions, 



B. Many particles 

After the one-particle eigenstates are obtained, these 
can be used as the single-particle basis for solving the 
eigenstates of the interacting many-body system by ex- 
act diagonalization. The A^-particle Hamiltonian can be 
written in the second quantization formalism as 



i,j,k,i 



(8) 



where Sj are the energy eigenvalues of the single-particle 
Hamiltonian, 



ijkl 



(9) 



are the matrix elements of the Coulomb two-body inter- 
action V{r,r') = l/\r — r'\ in the single-particle basis. 



and Vij = i^Tpi V{r) ipjj, where V{r) contains the Zee- 
man interaction and additional external potentials that 
are not included in Eq. (2). 

The interaction matrix elements of Eq. (9) can be cal- 
culated as follows. Let t/ji be the single-particle eigen- 
functions expanded in the sine basis of Eq. (7), 



V'i(r)=^<La(r) 



The interaction matrix elements are then 



1 



^Jkl = / dri / dr2**(ri)**(r2) — 4'fe(ri)*j(r2) 

R2 Jr2 ri2 

E °a*<"c"d / / dr2 



a.b.c.d 



X La{ri)Lb{r2) — Lc{r-i)Ld{r2) 
ri2 

a,b,c,d 



leads to 



ri2 



/ dri2 / d9 V (-^)"J„(fcrl2)e' 

Jo Jo „=_oo 

27r J dri2 Jo(A;ri2) = ^ . 



(fc) 



The potential l/ri2 can now be written as the inverse 
Fourier transform of 



1 

ri2 



as: 



ri2 



ri2 



(2' 



2tt 



— / dfc^d/c„-e'*^«(^^-^i)e''=''(''=-''i^ (10) 



With the identity of Eq. (10), the integrations over dif- 
ferent coordinates factorize in the interaction matrix el- 
ement: 



Vabcd= 



N 



27rL J^2 k 

oo 



/oo 
dxi sinc(xi — ax) sinc(xi — c^je^^^^^^ 
-oo 

/oo 
dyi sinc(yi - ay) sinc(j/i - Cy)e'''y'"' 
-OO 

X / dx2 sinc(a;2 — b^) sinc(x2 — dx)e~^^'^^^ 

J —oc 

X dy2 sinc(y2 - by) sinc(?/2 - rfy)e"*'=''^^(ll) 



The sine functions can be replaced by their integral rep- 
resentation 



1 

sinc(a;) = — / dt e'"^* , 
27r 
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and the integrals over x and y coordinates are of the form 

/oo 
da;sinc(a; — a) sinc(a:: — 6)e*'^^ 
-oo 

^e^'=»(27r - |fc|) , |A;| < 27r, a = 6 

, |fc| > 27r. 



By substituting this result into Eq. (11), the original four- 
dimensional integral over two planes reduces into a two- 
dimensional integral over a finite square in k-space, 



-2tv J-2tt 



r27r rK(e) 



= de dk 
Jo Jo 



X Ia.^c.^ {-k cos{0))laycy {-k sin(6')) 
xh^aAkcos{e))hydy{ksm{e)) , (12) 

where K{6) = 27r/ max(| 003(6*) |, | sin(0)|) is the radial 
integration limit corresponding to the square. The last 
form can be used in numerical calculations. 

One can see that in Eq. (12), one obtains five different 
integrals depending on how many of the four functions 
lab have the same indices. In addition, the case with 
two equal index pairs is naturally split into two cases, 
depending on whether the equal indices belong to the 
same Cartesian component of k. In most cases, some 
further analytic work can be done to handle the angu- 
lar integral. For instance, in the case when all the index 
pairs differ, such that ax ^ Cx, hx ^ dx, ay ^ Cy and 
by ^ dy, the integrand can be written as a sum of terms 
of the form cos {k[m cos{d) + nsm{9)]}, and the angular 
part can be integrated analytically, and we are left with 
a one- dimensional numerical integral. In this way, we 
are able to calculate the interaction matrix elements be- 
tween the sine basis functions, and then for any external 
confinement potential, Eq. (10) can be used to construct 
Vijki- 

It turns out that the calculation of Vijki from Eq. (12) 
is computationally very time-consuming, because one has 
to loop over foiir indices on both the right- and left-hand 
sides of Eq. (12). Luckily, this basis change can be triv- 
ially parallelized and a very efficient scheme can be ob- 
tained using graphics processing units (CPUs). 

We performed the calculation of the interaction ele- 
ments in Eq. (12) with an Nvidia Tesla C2070 graphics 
processing unit, which was programmed with CUDA^^, 
a parallel programming model for Nvidia GPUs. On the 
GPU, the computation is parallelized across tens of thou- 
sands of lightweight computational threads, which are 
organized in independent blocks. In our parallelization 
scheme, each block computes one element of Vijki ■ Inside 
the block, the sum over the index d is parallelized across 



the threads with each thread corresponding to a value of 
d. The threads then loop over the indices a, b and c, 
and in the end the results of all threads in the block are 
summed with a parallel prefix sum algorithm to obtain 
the final result. 

In Eq. (12), Wabcd does not depend on the state indices 
i, j, k, I, and it is beneficial to calculate it beforehand and 
store it in a table in the GPU memory. In double pre- 
cision floating point arithmetic, the size of the table for 
a, N X N mesh is 8A^^ bytes. For a 12 x 12 mesh, the 
size is approximately 3.3 gigabytes, which fits into the 
6 gigabyte global memory of the state of the art Tesla 
cards, such as the C2070. We also utilize the fast on-chip 
shared memory by caching the expansion coefficients a 
before the calculation. The GPU speeds up the matrix 
element calculation by a factor of around 13 when double 
precision arithmetic is used. 

Unfortunately, we have found the 12 x 12 mesh insuf- 
ficient for double quantum dot calculations if a realis- 
tic distance between the minima is used. Therefore, the 
calculation has to be divided so that the whole v^bcd 
matrix is not calculated at once. We lowered the mem- 
ory requirement by calculating first V^"^; by fixing the ax 
index in Eq. (12). As a consequence, it is sufficient to 
calculate the Uabcd also using a fixed ax index, and the 
memory requirement is dropped to 8iV^ bytes, which al- 
lows calculation with a 17 x 17 mesh. The Vijki elements 
are obtained by summing the V^"y elements over Gx- The 
sum is updated after the calculation of each V^j^^ matrix 
to save memory. This modification of the algorithm adds 
some serial work, which slows down the computation, but 
a compromise between memory requirement and speed 
must be made. 

The sum in Eq. (12) could be further divided to allow 
larger mesh sizes by fixing more indices, but the computa- 
tion time becomes fast a limiting factor. The calculation 
of interaction matrix elements for the 24 lowest single- 
particle states using a mesh size of 17 x 17 takes almost 
two days. The exact diagonalization part is much faster 
than this. 



IV. RESULTS 
A. Convergence of the single particle states 

In this section, we compute the single-particle eigen- 
states of systems consisting of 1, 2, and 4 minima using 
the Lagrange mesh method. The accuracy of the method 
and optimal parameters are discussed. In the following 
sections, we then use the obtained single particle states 
in the actual many-body computations. 

First, we studied the convergence of the method in the 
analytically solvable case of just one parabolic well. The 
confinement strength was fkvo = 4 meV. We computed 
the 24 first single particle energies with different mesh 
parameters TV and L and compared the results with the 
analytical Fock-Darwin eigenenergies. The relative dif- 
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FIG. 1. (Color online) The convergence of the 24 first single 
particle energies in the case of one parabolic dot as a function 
of A'' (the grid being A'' x A^). The relative differences of the 
single-particle energies (computed with the Lagrange mesh 
method) with the analytical Fock-Darwin energies are shown. 
The thin red curves show the first 24 states case with L = 280 
nm (the thick red curve shows the average relative difference). 
The thick dashed curve shows the average relative differences 
in the smaller grid case with L = 200 nm. 

ference of the energies can be seen as a function of the 
grid size N in Fig. 1. 

Fig. 1 shows that given large enough iV, the relative 
difference of the energies converges to the order of the 
numerical double precision accuracy in the L — 280 nm 
case. The effect of the size of the simulation area can also 
be seen in the figure. The smaller area case (the black 
dashed curve, L — 200 nm) shows faster convergence 
with respect to N . However, the finite simulation area 
results in some error as well, and thus the convergence in 
the L = 200 nm case stops before it reaches the double 
precision. 

The main topic of this paper is the simulation of 
singlet-triplet qubits. We will first study one-qubit dy- 
namics and then use our model to simulate a system of 
two singlet-triplet qubits. Next we discuss the conver- 
gence of the method in these systems. 

In the potential in Eq. (1), the derivative of the po- 
tential is not continuous; the min-function causes an edge 
at the interface of two branches. This sharp edge can be 
problematic in the Lagrange mesh method due to the fi- 
nite number of mesh points. To alleviate this, rounding 
of the edges was used in the case of multiple dots. The 
rounding was found to speed up the convergence of the 
single particle states. 

The rounding is achieved by defining a matrix R at 
each grid point. R has the different dot potentials in 
its diagonal. For example, in the case of four dots at 
locations ri...r4 the diagonal entries are Rn = Vi — 
^m*ujQ\r — Tip and so on. The non-diagonal entries are 
constant S and define the strength of the rounding. The 
potential at the particular grid point is given as the small- 



est eigenvalue of R. The effect of the rounding can be 
seen in Fig. 2. 
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FIG. 2. (Color online) The effect of the rounding on the DQD- 
potential. The potential is shown in the a;-axis. The minima 
are located ai x — ±40 nm. The confinement strength is 
fiojo ~ 4 meV. Both non-rounded (5 — 0, red dashed curve) 
and rounded (5 = 2 meV, black solid curve) potentials are 
shown. 

The current maximum grid size in the computation of 
the l^jfci -elements is = 17 due to the GPU memory 
limitations (larger grids can in principle be computed, 
but with the expense of considerably longer computations 
times). As the accuracy of the method depends non- 
trivially on both the simulation area L and the grid size 
N, the value of L was optimized. 

We compared the obtained eigenenergies with those of 
a large system (iV = 68 and L = 300 nm or L = 320 nm) 
and chose the value for L that gave the smallest error with 
respect to the more accurate large system. The relative 
difference of the energies as a function of L in the two 
dot case is shown in Fig. 3 and in Fig. 4 in the four-dot 
case. The potentials for the two- and four-dot systems 
are illustrated in the insets of Figs. 3 and 4. The two-dot 
potential consists of parabolic dots with the distance 80 
nm between their minima. In the four-dot system, dots 
1 and 2 are 80 nm apart, dots 2 and 3 are 120 nm apart, 
and 3 and 4 are 80 nm apart. 

Figs. 3 and 4 show that with = 17 the optimal 
value of L is between 180 nm and 200 nm in the two-dot 
case and between 260 nm and 290 nm in the four-minima 
case. Up to this point, the convergence of the energies is 
monotonous. With too small L, the wave function 'leaks' 
out of the simulation area, and with too high L the grid 
spacing becomes too large. The singularity like dips in 
the relative difference curves probably result from the 
fact that the errors due to finite L and N have different 
signs. At the dip, these errors nearly cancel each other 
out. 
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FIG. 3. The convergence of the 24 first Lagrange mesh single 
particle energies as a function of the simulation area length 
L in the two dot case. The two-dot potential is illustrated 
in the inset. The rounding is set to 5 = 2 meV. The grid 
size is 17 x 17 (A'' = 17). The relative differences of the 
Lagrange mesh energies with the energies of a larger system 
{L = 300 nm and A'' = 68, also computed with the Lagrange 
mesh method) are shown. The thick curve shows the average 
relative difference for the 24 states. 
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FIG. 4. The same as in Fig. 3 but for the four dot case, with 
the reference system being L = 320 nm and A'^ = 68. 

B. Singlet-triplet qubit 

In this section, we use the Lagrange mesh and ED 
methods to simulate the time evolution of the state of 
a singlet-triplet DQD qubit. We demonstrate that, by 
applying local electric and magnetic fields in our model, 
we can achieve full quantum control over the state of the 
qubit and reproduce realistic dynamics of the system in 
our simulation. 

We used a potential that consists of two parabolic dots, 

V{r) = iTO*womin{|r-ri|Mr-r2p}, (13) 



to model a singlet-triplet qubit. The confinement 
strength was fiwo = 4 meV and the distance between 
the dots was a = |ri — = 80 nm. Our DQD potential 
is illustrated in the inset of Fig. 3. 

The logical basis of a singlet-triplet qubit consists 
of the two lowest eigenstates, the singlet state, \S) = 
75 (I tl) - I It)), and the 5^ = triplet state, |To) = 
"^(1 ti) + I it)) (the arrows denote direction of the elec- 
tron spins). An arbitrary state of the qubit can be 
written as 

|V)=cos(^0|5)+sinQ)e^^|ro). (14) 

Here, G [0, tt] and </> e [0, 27r). The state of the qubit 
can thus be visualized using the surface of a Bloch sphere 
with \S) and |To) at the north and south poles, and 6 and 
<f) denoting the angles with respect to z and x-axes. 

It should be noted that the DQD is not a true two- 
level system. There are higher excited states as well, 
and Eq. (14) is just an approximation. However, in our 
simulations, the weighs of the higher states were found 
to be negligible with practical parameter values, and the 
system can be considered as two level in this sense. 

Universal quantum control of the qubit requires ro- 
tations around at least two different axes in the afore- 
mentioned Bloch sphere. In DQD singlet-triplet qubits, 
rotations around the z-axis are controlled by the ex- 
change interaction (the singlet-triplet energy difference) 
J — — Es and rotations around the x-axis can be 
generated by a magnetic field gradient ABz between the 
dots. The axis of the rotation in the Bloch sphere is then 

n = Jez+ ^BgAB^e^, (15) 

and the frequency is^ 

/ = ^J^ + ifiBgAB^y/h. (16) 

Here, g = —0.44 is the GaAs gyro magnetic ratio, /is the 
Bohr magneton and h the Planck's constant. 

The single particle eigenstates are computed using the 
Lagrange mesh method and they are then used in the two 
particle ED-calculations. In our model, the z-rotations 
are created by detuning (a potential energy difference e 
between the minima of the dots) the two parabolic dots, 
which lifts the degeneracy of the \S) and |To) states and 
results in exchange interaction. The x-rotations are cre- 
ated using a local magnetic field gradient that is taken 
into account by the Zeeman-term. 

The detuning potential and the local magnetic field 
are modeled as step functions that are zero far away 
from the dot minima and have different signs in the two 
dots. We calculate the the matrix elements Vij,cri,(j2 — 
(V'j.CTi |^|V'j.cr2), where V is cither the detuning po- 
tential or the Zeeman-interaction, in the eigenbasis 
{"01,(7} obtained using the Lagrange mesh method {a 
denotes the spin quantum number). The detuning 
and the Zeeman-term are then taken into account in 
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the two-body ED through the one-body operator V = 

The evolution of the initial state |V'(0)) of the qubit is 
computed by propagation, using 

mt + At)) = exp (^-^Hit)A?j (17) 

where H{t) is the (time-dependent) two-body Hamilto- 
nian. The matrix exponent is computed using Lanczos 
method. To study the evolution of the qubit 's state in 
the Bloch sphere, the angles 6 and (f) in (14) are extracted 
from \ip(t)) by using the properties of the two-body spin 
operator S^, i.e. S^\S) = and S^\To) = 2h^\To). 

The first 24 single-particle states were computed using 
the Lagrange mesh method. The mesh parameters were 
N = 17 and L = 210 nm. The rounding was set to 6 = 2 
meV. The Vijki- and V^ij-elements (corresponding to both 
the detuning and the Zeeman term) were computed for 
the 24 single-particle states. 

We first demonstrate the control of the qubit in a sim- 
ple case. In this simulation, the system is initially in the 
singlet state. The dots are detuned so that the differ- 
ence between their energy minima is e = V2 — Vi = 4.3 
meV. The detuning lifts the degeneracy of the singlet 
and triplet states, resulting in an exchange energy of 
J « 3.748 fieV. A magnetic field difference of AB2 = 0.4 
is then put between the dots and the system is let to 
evolve for 1 ns. The singlet and triplet probabilities were 
computed by projecting the state of the qubit onto the 
operator. 

The computed time evolution of the singlet probabil- 
ities p{\S)) can be seen in Fig. 5. Fig. 6 shows the 
evolution of the state of the qubit on the Bloch sphere. 

In Fig 5, the detuned singlet probability oscillates be- 
tween its maximum 1 and minimum 0.12. The singlet 
probability never goes to zero due to the z-rotation driven 
by the exchange energy J « 3.748 /leV. The frequency of 
the oscillation is / « 2.618 GHz, which is very close to 
the value given by Eq. (16), / « 2.625 GHz. In the non- 
detuned case, the probability oscillates between 1 and 0, 
as expected. In this case too, the computed frequency 
coincides very well with Eq. (16). 

Fig. 6 shows that in the detuned case, the plane of the 
rotation is tilted from the |S')-|To) plane, as expected by 
Eq. (15). The state never reaches \Tq) (the south pole) 
during the simulation. The non-detuned case oscillates 
between \S) and |To), passing through the spin localized 
states I ti) = l/V2{\S) + \To)) and | ;t) = i/V2{\S) - 
\To)). 

We also tried more complicated pulse sequences and 
tracked the evolution of the state in the Bloch sphere. 
One such is demonstrated in Fig. 7. Here, the detuning 
strength was oscillating, e{t) — eq sin(27r/t), where eg = 
4.6 meV and / = 1 GHz. A non-trivially time dependent 
J causes the axis of the states's rotation to change as a 
function of time, which leads to quite complicated paths 
on the Bloch sphere. 
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FIG. 5. (Color online) The time evolution of the singlet prob- 
ability p{\S)) (red curve). The detuning is e = 4.3 meV, and 
the magnetic field is AB^ — 0.4 T. This part is omitted from 
the figure, as p{\S)) is constant during it. The simulation 
time is 1 ns and the time step length At — 1 ps. The dashed 
cyan curve shows the non-detuned case with J ~ 0. 




FIG. 6. (Color online) The evolution of the state of the DQD 
qubit on the Bloch sphere. The qubit is initially in the singlet 
state (at the north pole). A detuning of e = 4.3 meV and 
a magnetic field gradient of AB^ = 0.4 are turned on and 
the state is let to evolve. The simulation time is 1 ns, and 
At=l ps. The thick red curve denotes the trajectory of the 
qubit's state on the sphere. The thin cyan curve shows the 
non-detuned case, with J « 0. The dashed black curve (the 
equator) shows the xy-plane, where p{\S)) = 1/2. 

In conclusion, our model, based on Lagrange mesh ED, 
can be used to simulate GaAs singlet-triplet DQD qubits. 
We can simulate the realistic full quantum control of the 
qubit starting from the first principles. Next, we proceed 
to use the model in studying two-qubit dynamics. 
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FIG. 7. (Color online) evolution of the state of the DQD 
qubit on the Bloch sphere. The qubit is initially in the 
singlet state (at the north pole). An oscillating detuning, 
e{t) = to sm{2n ft), where eo = 4.6 meV and / = 1 GHz, and 
magnetic field AB^ — 0.4 T, are then applied. The simulation 
time is f ns, and At = 1 ps. 

C. Two singlet-triplet qubits 

The Lagrange mesh allows the study of interplay of 
singlet-triplet qubits. For example, the entanglement 
of two singlet-triplet qubits by the long distance dipole- 
dipole interaction can be simulated using this method. In 
this section, we first compute and study the lowest eigen- 
states of the two-DQD system, using different detunings 
for the two qubits. The main topic of this section is the 
entanglement of singlet-triplet qubits. We show that our 
model can be used to simulate the entangling procedure 
demonstrated recently by Shulman et al^°. 

We model the two-DQD system with an external con- 
finement potential that is the minimum of four quadratic 
wells, 

V^(r) = imX'min{|r-r,f}. (18) 

Our simulation system can be divided to qubits A and B. 
A consists of the wells at ri and r2, with the dot distance 
o-A = \ — f2\ — 80 nm. Similarly, the inter dot distance 
of the qubit B is ub = {fs — = 80 nm. The inter 
qubit distance is given hy d = |t"2 — rsj — 120 nm. The 
confinement strength is hujo — 4 meV. The potential is 
illustrated in Fig. 8. The inter-qubit distance and the 
confinement strangth are large enough that there is no 
tunneling between A and B, so the qubits interact only 
through the Coulomb repulsion of their electrons. Also, 
the qubits interact mainly via the electrons in the dots 2 
and 3, as the inter-dot distances are quite large. 

The qubits A and B can become entangled due to the 
fact that under the exchange interaction, the charge den- 
sities of the \S) and \Tq) states differ. When the detuning 
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FIG. 8. Contour plot of the two-DQD potential. The system 
consists of four QDs, divided to two qubits A and B. The con- 
finement strength is huo ~ 4 meV. The qubit-qubit distance 
is d = 120 nm, and the distance of the dots in the qubits are 
tSA = as = 80 nm. The contours are shown in meV. 

lowers the potential energy in one of the dots of the qubit, 
the singlet state charge density becomes more located in 
this dot. However, if the detuning is not too high the 
triplet density is unaffected due to the repulsive exchange 
force in the spatially anti-symmetric triplet state. 

The singlet and triplet states have differing charge 
densities, and hence the Coulomb repulsion between the 
qubits depends on the states of the qubits. This condi- 
tioning creates an entangled state when the qubits are 
evolved under exchange. 

A bipartite state |'(/')ab & Ha <8) "Hb {Ha and Hb 
are the Hilbert spaces of the subsystems A and B) is 
an entangled state if it cannot be written as a tensor 
product \'4')ab = |V')a ® In general, if the vec- 

tor \iI})ab is written in any orthonormal product basis 
{\e^)A® \ej)B}ij, 

\i^)AB = ^M„-|e,)^(^ |e,)B, (19) 

it is an entangled state if and only if the matrix of coef- 
ficients, M = {My}, is not singular. 

The degree of entanglement can be determined by some 
entanglement measure. One such measure is the con- 
currence. In case of pure states, and two-level systems 
(qubits) A and B, concurrence C is given as C = 2VA1A2, 
where Ai and A2 are the eigenvalues of matrix M^M . It 
is easy to see that this simplifies to the formula 

C = 2|det(M)|. (20) 

Concurrence can also be generalized to mixed states^^. 

Concurrence assumes values between and 1. A non- 
zero C is a property of an entangled state, and the higher 
the value of C, the higher the degree of entanglement. 
The maximally entangled Bell states have C = 1. 
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In our two singlet-triplet qubit system, the Hilbert 
spaces are given as the two lowest eigenstates of 
a DQD-system, Ua ^ Ub = ITq)}. The 

M matrix is thus obtained by projecting the four- 
electron wave function onto the computational basis 
{|55),|5To),|To5),|ToTo)}; Mn = {SS\iP), M12 = 
(STol^), M21 = (To^lV) and M22 = {TqT^H). 

It should be noted that while the states \SS), \STo), 
and Ires') are eigenstates of the four-particle 5^ operator, 
\TqTq) is not. Indeed, it is not given as an eigenstate by 
the Lanczos iteration. In order to do the projections onto 
the computational basis, the state |ToTo) was generated 
using localized magnetic fields. 

As |To) = ^(1 n) - 1 m, \ToTo) = |To)a® |To)s can 
be written as 

\tqTo) = ^1 n)A 8> I n)s + ^1 n)A ® \ i-\)B 

+ \\i\)A®\n)B + \\i\)A®\i\)B. (21) 

In this decomposition, |ToTo) is written using the 5^ = 
eigenstates. These S^ = eigenstates can be generated 
using strong localized magnetic fields in the four dots 
of the two-qubit system. For example, | '\\) a ® \ ti) s 
is obtained as the ground state of a system where the 
magnetic field is up in the first dot, down in the second, 
up in the third and down in the fourth (the Zeeman- 
term alignes the spins of the electrons along the magnetic 
fields). 

Decompositions similar to Eq. (21) can be written for 
the other three states as well. Indeed, in the non-detuned 
case, the singlet states given by such decompositions were 
found to be the same eigenstates that Lanczos iteration 
would find. However, the aforementioned magnetic field 
scheme for creating the 5z = eigenstates can only be 
used to create states that have identical density in the two 
dots of the qubits. Hence, it is not well suited for creating 
the detuned singlet states. Fortunately if the detuning 
is in the practical operation regime of DQD-qubits, the 
|ToTo) density remains symmetric with respect to the two 
dots of the qubits. 

Thus, the computational basis can be created as fol- 
lows. The states \SS), |5To) and \TqS) arc given as eigen- 
states by Lanczos and they can be identified by their 
spin. The state |ToTo) is generated using the decomposi- 
tion Eq. (21). The wave function can then be projected 
onto this basis, and the concurrence can be computed 
according to Eq. (20). 

In the scheme where the qubits A and B are first 
brought to the xy-plane and then let to evolve under 
exchange (used for example by Shulman et al.^°), we can 
derive a simple analytic formula for the time dependence 
of the concurrence. 

In the absence of magnetic fields, the Hamiltonian of 
the two-qubit system is close to a diagonal one in the 
basis {\SS),\STq),\ToS),\ToTo)] (this was verified nu- 
merically). The diagonal entries of the projected Hamil- 
tonian are the energies Ess, Estq, EtqS and EtqTo- As 



the qubits are let to evolve in the xy-j>\anc, the weights in 
the M matrix obtain phase factors proportional to these 
energies. 

Let the system be initially in the state with Mij = ^ 
Vi,j e {1,2}, i.e. |?A) = I n)A ® I n)B- The system 
is then let to evolve. If we approximate the projected 
Hamiltonian to be diagonal, the time dependence of the 
coefficients My is given as Mij{t) — ie'^^j*/'', where 
Ell = Ess a-nd so on. Inserting these in Eq. (20) yields 
the formula for the concurrence, 

C{t) = ^^/2-2cos{AEt/h), (22) 

where, Ae = Ess + EtqTo - Estq - ExaS- 

In Eq. (22), the parameter A^; represents the cou- 
pling between the qubits. Eq. (22) shows that the 
entanglement indeed arises from the differences in the 
charge densities. If all the computational basis states 
have identical charge densities, the Coulomb repulsion 
between the two qubits, 7, is the same for all these 
states. In this case, the energies of the computational 
basis states are Ess = '^Es + 7, EtqTq = 2Eto + 7 and 
EsTo = EtoS = Es + Eto + 7. and = 0. The concur- 
rence is thus zero, i.e. there is no entanglement between 
the qubits. In the detuned case, the states have different 
densities and different values of the Coulomb repulsion. 
Hence, Ae j^O, and the concurrence oscillates according 
to Eq. (22). 

The single-particle states and the VijTj; -elements were 
again computed using the Lagrange mesh method. The 
simulation cell area was L = 280 nm, and the mesh size 
was N ~ 17. The rounding of the edges was set to S = 2 
meV. The 24 first one-particle states were used in the 
four-particle ED computations. As in the one qubit case, 
the detuning and local magnetic field matrix elements 
l^i,j,<Ti,CT2 = CTi l^lV'j.cra) were computed in order to 
achieve full quantum control over the qubits. 

First, we study the lowest eigenstates of the four- 
electron system. Without detuning (e^ = Vi — V2 = 
and €b = V4 — V3 = 0) the six first states (with Sz = 0) 
are close to each other in energy. 

The ground state is \SS) (s = 0), and the next two 
states arc |STo) and \ToS) (s = 1). The next three states 
given by Lanczos are what we call the triplet states, su- 
perpositions of s = 0, s = 1 and s = 2 eigenstates (in 
the four-particle case, there can be several S*^ eigenstates 
with given quantum numbers s and Sz) These three eigen- 
states of are so degenerate that Lanczos mixes them. 
The three triplet states share the same energy as |ToTo) 
(which also is not an eigenstate, but another linear 
combination of the triplet states). The electrons of the 
first six states are symmetrically located in the four dots, 
one electron in each. 

With non-zero detuning, one begins to see differences 
in the charge densities of the lowest eigenstates. Fig. 
9 shows the efi'ect of the sign of the detunings on the 
ground state \SS}. The charge densities of the lowest 
states (given by Lanczos) in the detuned case, ca = £b = 
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FIG. 9. The effect of the sign of detuning on the density of the ground state \SS). In the middle: the non-detuned \SS) density. 
The numbers in the middle plot refer to the dots 1-4 and A and B denote the qubits. The signs of the detunings, ea = V2 — Vi 
and es = V3 — V4, are shown by the axes and es. In the upper left corner: = 4.35 meV and en ~ —4.35 meV (dots 1 
and 3 have low potential). In the upper right corner: = es = 4.35 meV (dots 1 and 4 have low potential). In the lower left 
corner: ea = es = —4.35 meV (dots 2 and 3 have low potential). In the lower right corner: ea = —4.35 meV and es = 4.35 
meV (dots 2 and 4 have low potential). 



4.35 meV are shown in Fig. 10. 
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FIG. 10. The charge density of the lowest eigenstates of the 
DQD system. The detunings are ea = es = 4.3 meV (dots 1 
and 4 have low potential). The upper left plot shows the \SS) 
state, upper right the \STo) state, lower left the \ToS) state, 
and lower right the density of the triplet states (identical to 
the density of \ToTo)). 

There is a difference in the densities depending on 
which of the dots have the low potential, as can be seen 
in Fig. 9. In the upper right corner of Fig. 9 (and also in 



Fig. 10) the dots that are the furthest apart from each 
other have the lowest potential. This facilitates the lo- 
calization of the singlet state in these dots, as it reduces 
the Coulomb repulsion. In the cases when dots 2 and 4 
or 1 and 3 are in the low potential, the singlet localizes 
only in the further away dots, 1 and 4. In the case where 
dots 2 and 3 have the low detuning, there are four identi- 
cal peaks, the singlets cannot localize in the neighboring 
dots due to the Coulomb repulsion. 

The dipole-dipole entanglement effect relies on the dif- 
ferences of the singlet and triplet densities. Thus, to 
this end, the optimal detuning configuration should be 
as in Fig. 10, the dots furthest away are detuned to low 
potential energy. The densities in Fig. 10 show that 
the singlets localize to the dots that have lower poten- 
tial (dots 1 and 4). The fourth plot represents all the 
triplet states and is identical also to the \TqTo) density. 
It shows four identical peaks, with exactly one electron 
in each dot. When the detuning is further increased, the 
singles localize fully to the low lying dots, and at very 
high detunings, the triples start to localize as well. 

With high detuning (e^ and cb above 5 meV), \SS), 
jS'To), and \TqS) were still lowest in energy. However, 
the triplet states were not the next three in this case. 
There are states lower in energy than the triplets, in- 
cluding other instances of the states 155), \STo) and 
|ro5) (i.e. states that are of the form \X)a €5 \Y)b, 
where \X) and \Y) are s = or s = 1 eigenstates of 
the two-electron 5^-operator). Next we consider states 
|55), \STo), |To5), and \ToTo) that have the lowest en- 
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FIG. 11. (Color online) |A_e| as a function of the detunings e,4 
and e_B- The values of |A_b| are shown in meV. A high value 
of |A_e| means fast qubit-qubit coupling. The detunings are 
sampled at intervals of 0.2 meV. 



ergy and study the dependence of the coupling parameter 
A_B = Ess + EtoTo - EsTo - Et„s on the detunings 
and es was studied. A high value of \Ae\ means fast 
qubit-qubit coupling, as seen in Eq. (22). In Fig. 11, 
jA^;! is shown as a function of the detunings. 

The area of high lA^j {eA,£B between 4.3 meV and 
6.5 meV) is roughly rectangular. Because the transitions 
from 5'(1,1) to 5(2,0) and from ro(l,l) to ro(2,0) hap- 
pen quite abruptly with respect to increasing detuning. 
In the low detuning region, the singlet and triplet states 
are both localized in two dots of the qubits and in the 
very high detuning case even the triplets localize to dots 1 
and 4. In the rectangular high jA^jj-area, the singlets are 
in (2, 0)-configuration and triplets in (1, 1) configuration. 

Outside of the main rectangular area, there are also 
areas of smaller increase in | A^; | . These are located at the 
sides of the rectangular peak, and probably result from 
the fact that one of the qubits is in the detuning interval 
4.3 meV and 6.5 meV. These side areas were also present 
in the computations done by Nielsen et al, in addition to 
a large plateau in |A£;| when both detunings are high^^. 
This plateau is not present in our results, possibly due 
to the fact that our qubit-qubit distance is quite large. 

We then use our singlet-triplet qubit model, introduced 
in the previous section, to simulate the entanglement of 
the qubits A and B by the dipole-dipole interaction. We 
start from the ground state of the system, \SS) with zero 
detuning. The magnetic field gradients in the qubits are 
then turned on, and the state of the system is let to 
evolve. When the qubit reach the xy-j>\ane in the Bloch 
sphere, the magnetic field is turned off, and the detunings 
are turned on. When the detunings have reached their 
maximum values, the system was is to evolve again. The 
evolution is computed according to Eq. (17), and the 
concurrence according to Eq. (20). 



The computed concurrences can be seen in Fig. 12. 
Here, the qubits are first brought to the xy-plane using 
magnetic field gradients of AB^^a — AB^^b — 10 mT. 
We study the effect of the speed by which the detunings 
are increased to the values = es — 4.3 meV. Cases 
of an instantaneous increase and an adiabatic increase 
during a time of 1 ns can be seen in Fig. 12. The qubits 
are let to evolve in the xy plane for a time of 10 ns, and 
the concurrence is computed at each time step according 
to Eq. (20) (the time step length was At = 1 ps). Fig. 
12 also shows the concurrence given by formula Eq. (22), 
where the energies are computed by Lanczos (-EtqTq was 
obtained as the eigenenergy of one of the triplet states) . 
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FIG. 12. (Color online) Entanglement of two singlet-triplet 
qubits measured by concurrence. The qubits were brought 
to the a;j/-plane of the Bloch sphere and let to evolve for 10 
ns. The detunings were tA = ts = 4.3 meV. The figure shows 
only the evolution in the a;j/-plane, as the concurrence was zero 
before this. Solid red line shows the curve that was obtained 
with adiabatic increase of the detunings during 1 ns. In the 
dashed curve, the detunings were turned on instantaneously. 
The curves are obtained by projecting the wave function onto 
the computational basis. The black markers show the con- 
currence given by the analytical formula Eq. (22). 

In the adiabatic case, the concurrence oscillates be- 
tween and 1 as the state of the qubits is evolved in 
the xy-plane. The obtained curve coincides almost com- 
pletely with the one given by Eq. (22). In the non- 
adiabatic case, the frequency of the concurrence oscilla- 
tions, given by A^;, is the same as in the adiabatic one, 
but the amplitude is smaller, only about 0.81. The reason 
for this is that when the detunings are increased instanta- 
neously, the wave function leaks from the computational 
basis. Indeed, the the probability for being in the com- 
putational basis, Yld j=i l-^y Pi is only 0.82 in this case 
(during the evolution in the a;y-plane). 

The forms of the computed concurrence curves are sim- 
ilar to the measurement by Shulman et al.^^ apart from 
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the dccohcrcncc effects that have not yet been imple- 
mented in our model. Fig. 12 shows that the formula in 
Eq. 22 indeed describes well the entanglement of qubits 
in our system if the detunings are increased gradually. It 
is however limited to the case of a;?/-rotations. In the more 
complex cases of time-dependent detunings and magnetic 
fields, the concurrence can be computed by projecting the 
wave function onto the computational basis. In order to 
obtain the maximal degree of entanglement, the detun- 
ings should be increased adiabatically to their maximal 
values. 

It should be noted that it is not certain that the 24 
first single particle states are enough to produce quan- 
titatively accurate results in this two-DQD case when 
the detunings are very high. We have not studied the 
convergence of, for example, the lowest four-particle en- 
ergies by comparing the 24-basis results to ones obtained 
with a larger basis. Hence, the analysis in this section 
should be considered only qualitatively valid. One could 
of course compute more states, but increasing the single 
particle basis size will result in significant increase in the 
four-particle basis used in the computation of dynamics. 
The full dynamics computations were quite slow in the 
24 state case already, and as the main topic of this article 
is the presentation of our method, quantitative analyses 
are not included. 

Instead of doing the matrix exponent for the full four- 
particle Hamiltonian, one could first project it to some 
number of eigenstates of the Hamiltonian with zero de- 
tuning and magnetic field gradients. This small matrix 
can then be diagonalized exactly. This method was tested 
in the case of the 12 first four-particle eigenstates (the 
number of single particle states in the Lanczos was 24), 
and it gave the same results as the full dynamics com- 
putations. By using this projection method, one could 
increase the single particle basis size to obtain more ac- 
curate results without the expense in simulation lengths. 



V. CONCLUSIONS 

In summary, we presented a Lagrange mesh based 
scheme for studying many-particle states in lateral quan- 
tum dots. We introduced a Lagrange mesh based method 
for many-body systems and then proceeded to use the 
method in simulating singlet-triplet qubits. We intro- 
duced a model for simulating the full quantum control of 
singlet-triplet DQD qubits, and showed that this model 
can be used to produce realistic dynamics of the qubit 
system. The singlet-triplet qubit model was then used to 
study a system of two qubits. We computed and studied 
the lowest eigenstates of this system and also discussed 
the effect of electrostatic detuning on the eigenstates. 
The entanglement of the two qubits via the dipole-dipole 
interaction was simulated. 

The Lagrange mesh provides a very flexible method 
of dealing with complex confinement potentials in ED 
calculations. It allowed us to create a realistic first prin- 



ciples model for studying the interplay and dynamics of 
two singlet-triplet qubits. This method could be used to 
study complex effects that are difficult to include; in sim- 
pler models. In addition, our model could quite easily 
be further improved by the inclusion of the decoherence 
effects due to the environment, like in^^'^^. 
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APPENDIX: SINC MESH AND GAUSS 
QUADRATURE 

In this appendix, we consider the validity of the Gauss 
quadrature approximation of the potential energy matrix 
elements. 



Va'a « V{a)Sa'a, 



(23) 



in the case of the sine mesh. We mainly concentrate on 
the case of a parabolic potential, as this paper deals with 
locally parabolic quantum dots. General results concern- 
ing other potentials are briefly discussed. 

In principle, the potential matrix elements Va',a can 
be computed analytically. 



Va'. 



f 

J — ( 



ArLa'{r)V{r)La{r). 



(24) 



The approximation (23) can be considered valid if the 
results (i.e eigenstates and energies) obtained with it do 
not deviate considerably from the ones obtained with the 
analytical formula. 

Consider a one dimensional sine mesh over the interval 
(— -J, -j). The number of basis functions is N and the grid 
spacing is /i = ^. Unfortunately, the analytic integral 
for the matrix elements, 



1 r 



sine 



Xa') V{x) sine —{x — Xa 
. (25) 

is divergent for example in the case of a parabolic po- 



tential, V{x) 



In Eq. (25), 



+ ha and 



Xa' = + ha'. The integration can, however, be done 
over some large interval (— M, M), where M is a multi- 
ple of h, M = Nh and N > N. In this case, with the 
parabolic potential, the elements are given as 



Va'a =xlda'a + i-l 



,a +a 



2hM 



-M+Xa 



+ 



M-Xa 



t 



dt, (26) 
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where 5a'a is the Kronecker delta. The third term of (26) 
can be approximated as 



-M+Xa 



M- 



sm 



di 



< h min log 



N 



(27) 



oo or /i — )■ 0. Unfortunately, 
2hM „ 2L TV generally diverges, 



TT 



This vanishes when N 
the second term of (26), 
as N approaches infinity (i.e. we approach the analytical 
formula (25)). However, this also can converge to zero in 
the hmit A'^ oo if ft, at the same time approaches zero. 
This is the case for example if iV = ceihng (logiV) TV. 
Thus, given small enough h, the potential energy matrix 
obtained analytically is arbitrarily close to the diagonal 
one obtained with the Gauss quadrature. 

In the two dimensional case, the analytical integrals 
reduce to the same one dimensional integrals discussed 
above. However, due to computational limitations N 
cannot be very large in 2D when computing the two body 
interaction matrix elements (iV = 17 is the current max- 
imum for the computation of the Vijfci-elements because 
of the GPU memory limitations) . With realistic values of 
L, the off-diagonal terms of (26) are non-negligible with 

N = n. 

Consider for example the case of a parabolic well, 
V{r) = ^m*uiQr^. The values of the constants are 
m* = 0.067me, huo = 4: meV, L = 200 nm, L = 6L 
and N = 17. In atomic units, the absolute value of the 
off-diagonal elements is approximately 0.5 (the maximum 
being 0.501), while the maximum diagonal value is 7.184 
(The max;imum diagonal value contains the maximum 
value of the potential, 



,2 r2 



L . Most of the diagonal 
vahies are much smaller than this). Thus, the diagonal 
Gauss quadrature approximation does not seem to hold. 

However, the eigcnencrgies computed with this ana- 
lytically obtained potential matrix arc very close to ones 
obtained with the approximation. For example, the rel- 
ative differences between the 24 lowest eigenenergies are 
all below 10^'' (the energies from the approximation are 
also very close to the exact solution of the problem, the 
relative differences being again less than 10~^). Similar 
accuracy holds for the eigenfunctions. Higher the energy 
the bigger the differences are, but nevertheless the ac- 
curacy of the Gauss quadrature seems quite unexpected. 
Next, we try to find reason behind this phenomenon. 

For simplicity, consider again the one dimensional case. 
Similar arguments apply to the 2D case. Let u be an 
eigenvector of the Hamiltonian H, where the potential 
matrix is computed with the Gauss quadrature approx- 
imation, Hu = Eu. In order for the approximation to 
be valid, u should also be an approximate eigenvector of 
the Hamiltonian H, where the potential energy matrix 
is computed analytically, according to (26). Here, we set 
the interval of integration L to be large enough that the 



term in (27) can be neglected, N » N. Now, H can be 
written as 



H = H + A,A„.„ = (-1) 



a+a 



(28) 



where 5 = 

The symmetric matrix A is of rank 1, its only non-zero 
eigenvalue being = NS. The corresponding eigenvec- 



tor IS 
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-1) 



JVl 



— _L X ... v^— J - By projecting u onto the 
(mutually orthogonal) eigenvectors of A, one can com- 
pute the product Au. As, apart from A^, all eigenvalues 
of A are zero, it holds that 



flu = Eu + Au = Eu + (5(v'^u)v, 



(29) 



where v = [l — 1 1 ... (—1)-'^], and v'^u = ui — U2 + 
U3- Ui + ... + {—1)'^un. 

We now reason why v^u should be small. By 
the Lagrange condition, Ua = htj){xa), where ip{x) = 
UaLa{x). As i/i is a finite sum of diffcrentiable func- 
tions, it also is a diffcrentiable function. Thus, u is a 
discretization of a diffcrentiable function. Hence, the in- 
crements Ua —Ua+i should be small and cancel each other 
out, as (with large enough L) ui « and un ~ 0. 

This argument can be further justified by writing Ua = 
u{xa), and noticing that v'^u is related to the integral of 
the derivative of u. 

T _ ujxi) - u{xi + h) ^ ^ u{x3) - u{x3 + h) ^ ^ 
V u — - ii -f- - n ~\~ ... 

h h 

« --Y,u'{x2j+i)2h 



\ J"^^'^ u'{t)dt = \ H-L/2) - u{L/2)] . 



If the system area L is large enough, ip{±L/2) 0. Fur- 
thermore, as u{x) = htp{x), u{—L/2) — u{L/2) can be 
considered to be very close to zero. 

Now, we can approximate the error term Au using 
an error formula for Riemann-sums. As we have ap- 
proximated the ratio of differences, ^ of u{x), to be 
the derivative of m, some error arises from this approx- 
imation as well. From the Taylor expansion: — 
u\x) + 0{h^) = htfj'{x) + 0{h^), where we have used 
the fact that u{x) = hijj{x). The error for the Riemann 
sum formula is in this case smaller than -^p-, where 
D = maxx \tp'{x)\. The elements of the 'error' vector 
Au = 5(v^u)v thus obey 

\{^uU<^-^+SO{h') = ^ + 0{N-% (30) 



7V2 



27V4 



_We have shown that the elements of Au behave as ~ 

. They become small even with relatively small values 
of N. Eq. (30) also explains why higher eigenenergies 
tend to differ more between the Gauss approximation 
and the analytical formulas. Higher eigenstates oscillate 
more and thus have larger values of |V''(2^)I- 
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Same kind of analysis can be done to potentials V(x) = 
x'^,k G N OY V{r) = r'^, fc e N, and the results are simi- 
lar; the analytical formula for the potential energy matrix 
elements gives approximately the same eigenstates as the 
Gauss quadrature. The matrix A is of the same alter- 
nating form in all of these cases as in Eq. (28), and thus 
the same convergence arguments apply here. In princi- 
ple, the result can then be extended to any potential that 
can be written as a power series. 



In conclusion, the analytic formula for the potential 
matrix elements of a quadratic well can be made to con- 
verge to the Gauss quadrature approximation when the 
grid spacing goes to zero. In addition, the accuracy of the 
Gauss quadrature eigenenergies and states remains good 
even when the approximation does not hold for the indi- 
vidual matrix elements. Similar results apply for other 
potentials as well, for example r*^, fc e N. 
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